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Abstract 

The algebraic irrational number golden ratio (p = (1 + V5)/2 = one of the two roots of the algebraic 
equation x 2 -x-l=0 and the transcendental number n = 2 sin -1 (1) = the ratio of the 
circumference and the diameter of any circle both have infinite number of digits with no apparent 
pattern. We discuss here the relative merits of these numbers as possible random sequence sources. 
The quality of these sequences is not judged directly based on the outcome of all known tests for the 
randomness of a sequence. Instead, it is determined implicitly by the accuracy of the Monte Carlo 
integration in a statistical sense. Since our main motive of using a random sequence is to solve real 
world problems, it is more desirable if we compare the quality of the sequences based on their 
performances for these problems in terms of quality/accuracy of the output. We also compare these 
sources against those generated by a popular pseudo-random generator, viz., the Matlab rand and the 
quasi-random generator halton both in terms of error and time complexity. Our study demonstrates 
that consecutive blocks of digits of each of these numbers produce a good random sequence source. 
It is observed that randomly chosen blocks of digits do not have any remarkable advantage over 
consecutive blocks for the accuracy of the Monte Carlo integration. Also, it reveals that n is a better 
source of a random sequence than q> when the accuracy of the integration is concerned. 

1. Introduction 

The distribution of digits in an irrational number 1 — both algebraic and transcendental — has not 
been extensively explored. Only to a limited extent this has been discussed in [1], So far as a rational 
number is concerned, it could have either infinite number of digits in the conventional decimal 
number system or it could have a finite number of digits in a number system (not necessarily 
decimal). Whenever a rational number has infinite number of digits, these digits are necessarily 
recursive in some form(s). Hence the distribution of its digits is usually known/can be determined or 
even predicted. Our conjecture is that digits of a rational number cannot be uniformly randomly 
distributed. Much effort is not required to prove/disprove this conjecture. However, our 
inquisitiveness is with respect to the distribution of digits in an irrational number in which the 
knowledge of the current and previous digits does not help to predict the next digit. 
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1 An irrational number is one which cannot be expressed as the ratio of two integers. This implies that an irrational 

number will have infinite number of digits for its exact representation. 



To start with we have chosen one popular irrational number, viz., golden ratio q> from among the 
algebraic numbers and one extremely used more popular irrational number, viz., n from among the 
transcendental numbers 3 . We have studied the distribution of the digits of <p [1] as well as those of 
n , which are then employed as a random sequence source (RSS) to compute the values of Monte 
Carlo integrals. The question “Which is better — (p or n — in uniformity in distribution of its 
digits as well as in computing an integral?” is attempted for an answer. In fact, similar exploration 
could be carried out for other well-known/famous irrational numbers or even ordinary irrational 
numbers. Some of the ordinary irrational numbers could possibly turn out to be 
remarkable/outstanding in terms of uniformity of their digits or randomness of their digits. Some of 

the other famous irrational numbers are the Hilbert number 2^ « 2.66514414269023 , the Euler- 

Mascheroni constant y = lim — In «) ® 0.57721566490153 , and the numbers 

*-i k 

i‘ =e-’ 1 ' 2 *0.207879576350762, n e « 22.459157718361 1 (believed (not proved) to be a 
transcendental number) and e” « 23. 1406926327793. Such an exploration will help one to decide 
which irrational number should be used as an RSS not only for Monte Carlo integration but also for 
randomized/evolutionary algorithms such as ant system approaches, genetic algorithms, and 
simulated annealing. Some randomized algorithms perform better with a quasi-random sequence 
while others perform better with a pseudo-random sequence [2-36]. It may be noted that quasi- 
random sequences are more uniformly distributed (with less discrepancy) than pseudo-random 
sequences. One important advantage of using an irrational number for an RSS is that it obviates the 
need of employing a random number generator. We would need, for a reasonable real world 
problem, a generation of truly large random sequence — a sequence of billions or possibly trillions 
of random numbers — to solve multidimensional problems. Here the extra time of generation of 
random numbers are eliminated since we can simply pick out blocks of digits consecutively from an 
irrational number. However, such an extra time appears near about the same as that required for 
retrieval of numbers. An efficient storage and retrieval system possibly could improve the time 
complexity and could even prove better than instant generation of random numbers. It may be 
pointed out that not all irrational numbers are having digits randomly occurring. For instance, the 
Liouville number® 0.1 100010000000000000000010000 which has a 1 in the 1 st , 2 nd , 6 th , 24 th , 120 th 
etc. places and 0s elsewhere, is an irrational number where 0’s and l’s are not random. These are 
exactly predictable. 

In section 2, we describe simple Monte Carlo procedures for numerical integrations by 
appropriately slicing the n-dimensional domain after mapping it onto a standard domain. Since our 
main purpose of using a random sequence is to solve a real world problem, it is more desirable if we 
compare the popular random number generators along with the concerned irrational numbers based 
on their performances in terms of quality (error bounds) and cost (computational/time complexity) of 
the solutions that they produce. In section 3, we present the numerical results considering typical 
single integrals. Section 4 comprises conclusions. 


2 An algebraic number is a root of any rational polynomial. 

3 Transcendental number is a number that is not the root of any polynomial (of finite degree) whose coefficients are 
rational (or integer) numbers. That is, it is not the root of any integer polynomial implying that it is not an algebraic 
number of any degree. Every real transcendental number must be irrational since a rational number is, by definition, an 

algebraic number of degree 1. Some of the proven transcendental numbers are In 2, e, Jt, 2^, sin(l), r(l/3), and 



2. Monte Carlo Integration 


A Monte Carlo (MC) procedure is a method based on using uniformly distributed random numbers 
(RNs). This procedure is sometimes preferred over a numerical quadrature when the integrand is 
violently fluctuating or involves long trigonometric/special functions. We have chosen here the 
Monte Carlo integration for the purpose of comparing n and cp as RSS’s. Also, compared are n 
and (p against the popular pseudo-random generator Matlab rand and quasi-random generator 
halton. 

Let the single integral be 

b p 

I = \f(x)dx = \f(x)dx , (1) 

a a 


where the function f(x ) is continuous, non-negative, and single real-valued in the interval [a, V\. 
The MC procedure is as follows. 



Figure 1 MC integration of f{x) in [a, b] determines the area under the curve y - f(x) from 
x = a to x = b . 

S.l Choose M > largest value of f(x ) in the interval [ a,b ] . Initialize hit = 0 and N = 5000 (say). 

S. 2 Generate/choose a pair of random numbers (x,y) such that xe[a,b] and y e [0,M] . If 
f(x)< y then it is a hit and hence update hit = hit + 1 . 

S. 3. Continue S.l N times. The integration value is then I = - a )M ■ 

This MC procedure can be generalized for the integrand / (*) which has both negative as well as 
positive values in the interval [ a,b ] by introducing the minimum value R of f(x ) , replacing 
y e [0,M] by y e [ R,M ], and appropriately determining hit with proper signs. An alternative way 
is to divide the integration into or more parts such that in each part the function f{x) is either 
wholly nonnegative or wholly negative. WE then apply for each part the procedure similar to the 
foregoing MC procedure. 

However, we prefer to use MC method on an integral with fixed limits of integration, say, 0 and 1 
since the summation-based MC method (presented later) remains invariant. Any integral with limits 
of integration [a, b] can be mapped onto the integral with limits of integration [0,1] where the new 
integrand is a linear transformation of the original integrand. For example, the integral 



I = ]f{x)dx=\F{t)dt. 

a 0 

where the new integrand F{t) = pf (pt + q), p = b- a, q = a, dx = pdt . These values are derived 
from the linear transformation x = pt + q . When x = a,t = 0^>q = a. When 
x-b, t = l=> p — b — a .If o = 0, 6 = oo, then 


oo r oo 

I = | f(x)dx= jf(x)dx + \f(x)dx , 

0 0 r 

where r is a finite positive real number. The interval [0, r] can be changed, as before, to [0,1] . For 
the interval [r,oo], use the nonlinear transformation x-\ /(pt + q). The interval [- 00 , 00 ] can be 
similarly changed (mapped) to finite intervals [37]. 

Without any loss of generality, we now consider the single, double, and triple integrals 

A = \f(x)cbc, I 2 = jj ' f(x,y)dydx, / 3 = j j jf(x,y,z)dzdydx . 

0 00 000 

The unbiased MC estimate of the values of /, , I 2 , / 3 are 


7 i =-£/(*, )> 7 2 ^“jZZZ f( x »yj’ z k), 


r Ti 


,=i ,= 


=1 j = 1 *=1 


where each of x j , y p z k is chosen uniformly randomly in the interval [0,1] . The constant r = N in 
the expression of /, is the number of random numbers used to compute /, . The constant r 2 - N in 
the expression of / 2 is the number of random numbers used to compute I 2 while the constant 
r 3 = N in the expression of / 3 is the number of random numbers used to compute / 3 .The 
generalization to an n -dimensional integral is straight-forward. The*MC estimate tends to the true 
value of the integral as N —> 00 . We, however, restrict ourselves to single integrals for the purpose 
of comparison between n and (p . The MC multiple integration will be discussed in a future work. 

Let the function f(x) be square integrable 4 . From the central limit theorem, the standard error 

N 

Z(/0O-/,) 2 > a °d an estimate 

1=1 

of the error in /, is E = Js / N . However, a general numerical way of the relative error estimation 
[38] is to compute /, with N = 5000 uniformly distributed RNs (pseudo-RNs or preferably quasi- 
RNs) assuming sufficiently large precision RNs. Call the resulting integration value as the quantity 
of lower order accuracy Q' . Once again compute /, with N = 50000 uniformly distributed RNs. 
Call the resulting integration value as the quantity of higher order accuracy Q. The relative error is 


tends to XI 4n , the variance of the MC estimate I, is S = 

1 N - 1 


00 

4 A function is square integrable if || f(x) \ 2 dx is finite. 

—00 



then (Q - Q')/Q ■ If the precision (length) of the RNs is not sufficiently large then even with higher 
number of RNs, accuracy may not improve. It may sometimes even become worse due to clustering. 

Dividing the interval/domain of integration For the sake of obtaining desirable accuracy in the 
integration value, it is, in general, necessary to keep a provision of writing the given integral into a 
series of integrals by appropriately slicing/dividing the interval/domain into a number of 
subintervals/subdomains. Integrate the function in each subdomain and add the resulting integration 
values. The subdomains may or may not be equal. We, however, have considered them equal 
irrespective of whether the function (integrand) is ill-posed or not in a subdomain. For a real-world 
problem, such a programming simplification, in a Matlab environment or for that matter in most 
other programming environment, does not have a significant negative effect on the computational 
complexity or on accuracy. Since the sensitivity of the integrand could vary from one subdomain to 
another, accuracy of computation will also vary from one subdomain to another. When we add up 
the resulting integration values with unequal accuracies, the final integration value will be less 
accurate than the least accurate integration value. However, when one works with a precision of 15 
digits with appropriate equal subdomains, the relative error propagation is most unlikely to hit the 
fifth digit in the final integration value. It may be noted that a real world engineering implementation 
cannot be usually more accurate than four significant digits. This is because no measuring device 
can measure, in general, a quantity with more than 0.005% accuracy [38]. We can write, allowing 
a = a x and p = a n , />,. = a,. - o,._, , q t = a,_, , a , > a M , i = 2(1 )n - 1, (i.e., i = 2,3,4, •••,«- 1, ) 

I = \f(x)dx= J f(x)dx = ^ \f{p,t + qfp.dt 

a a } ,= 2 o 


If/(x) = xsin(l/x), a = a, and P = a n , p t =a t -a,_ x ,q i = a w , a, / = 2(1)«-1, then 

P n-\ 1 J 

/ = \x sin(l / x)dx = £ j(a, - )[(«, - a iA )t + a t _, ] sin(- )dt . 


If f(x ) = S>n - , a = a, =0 and P = a n = oo , then, using the transformation x = 1 /(pt + q) [37], we 
x 

have 



- — sin (- — )]dt = — (exact value). 


Considering the second integral, we have a = 0, P = 1 if we desire to change the interval of 
integration [0,oo) to [0,1]. However, we can instead choose oo as 20000 (say) for sufficient 
numerical accuracy for real world implementation and write 


20000 


l^dx, \[^dx = 20000\- 


sin(20000/) , n 
■at ~ — . 


20000 / 2 

u u u 

It may be observed that as t —> 0, sin(20000/)/(20000/) — > 1 (in the limit). Or, in other words, if we 
choose / = 10 -10 , we automatically get, in Matlab, sin( 20000/) /(20000/) as 0.999999999999333 



which is correct up to 12 decimal digits (sufficiently accurate for a real world application). Thus we 
may replace the lower limit of integration, viz., 0 by 10" 10 whenever a division by zero could occur. 
We may also vary the decimal number 20000 and observe the change (if any) in the value of 
integration up to, say, seven digits. 

3. Numerical Experiments 

We have considered several typical single test integrals having variable sensitivity over the interval 
of integration. Also, we have considered an integral whose solution is not known. Our focus is to 
compare n and <p to determine which one is a better source of a random sequence although we have 
brought in the pseudo-random generator, viz. the widely used Matlab rand generator and the quasi- 
random generator, viz., the halton generator [2, 3, 1 1]. The quasi-random generators for Monte Carlo 
integration are expected to perform better since these produce more uniformly distributed (with less 
discrepancy) random sequences than those produced by pseudo-random generators. Our numerical 
experiments also depict this fact. 

Let the general form of our single integral be as given in the equation (1). We have used the 
following notations in Table 1 that depicts Monte Carlo integration values for the following single 
integrals // and I 2 along with relative error in each of rand, halton, pi, and phi: RSS = random 
sequence source/generator, n = number of subintervals for the interval (/?-«), N = 5000 = 
number of random numbers used in each subinterval, True value = True/exact value of the integral. 

The single integrals along with their true/exact values, that we have considered are 


(i) /, = J x 2 dx = )-, 

0 ^ 

OO • 1 , , 10000 • 

_ fSinx , r r sin:t 1 . / 1 ... , n r sinx , 

(ii) I 2 = I dx= |[ + - sin(- )]dx = - * I dx, 

o x J X l-x l-x 2 l 0 i l0 x 

1 03 • 10000 

(iii) 7j = J(^^) 2 dk = — » | (^^) 2 afr S 1.5707963267949 , 

o * ^ io -10 x 

03 4 I j i 

(iv) I A = f jdx= f[ =- + ]cb C = -*l.. 

4 h + x 2 1 + JC 2 (1-Jt) 2 +1 2 


5707963267949 , 


10 " 


VI VI 

(v) / 5 = Ltsin— dx& xsin—dx « 0.3785283345 

J ^ J x 


(the value 0.3785283345 of the integral I 5 is computed dividing the interval [10 , 1] into one 

million equal subintervals since the exact/true value of I 5 is not known.), 


10000 


T rCOSJx , _ rCOSiJC , 7t „ ”f’COs3x, 

(vi) L - I -- - "--d x - 2 dx ~ — T « 0.00389361481 w 2 , d x , 

i* 2 +4 J jc 2 +4 2e 6 l jc 2 +4 

03 03 10000 

r fjtsinjt , * rjcsinx , n , , ^ r xsmx , 

(vii) L = T-dx = 2 -dx = — * 1.155734979 *2 -rdx, 

il + jc 2 0 J l + x 2 e 0 J l + x 2 

(viii) / s = jlnxdx = [jclnx- x], -[jclnjc-jc] 10 . 1(( « -.999999997597415 , 


10 “' 



(ix) g(a>) = 


I 


cos me sin* 

x 


dx = \ 


n H,Q < a> <\ 
n!A, co = 1 
0 , w> 1 


We take o> = 0.5, 1,2 for the comparison. Our 7, = g(0.5), 7 10 = g(l), /,, =g(2). Also as in the 

earlier integrals, we choose the limits of integration as [10' 10 , 10000] for the actual limits [0, qo] . 

Tables 1 and 2 provide us the MC integration (MCI) value of each of the integrals when pseudo- 
random generator Matlab rand, quasi-random generator halton, RSS’s pi {n) and golden ratio {(p) 
are used. Wherever the limit of integration zero causes division by zero in a numerical computation, 
zero(0) is taken as 10' 10 , otherwise zero is kept as such. For the limit of integration oo, we have 
taken 10000. This is expected to provide four significant digit accuracy which is, in general, 
acceptable in most applications. It may be noted in this connection that a measuring device cannot 
measure usually an accuracy greater than 0.005% [38]. However, for an intermediate integration 
value used for further computation, we should compute the value with higher accuracy/precision so 
that the final output to be used in real world implementation has a 0.005% accuracy. 

Table 1 Monte Carlo integration (MCI) along with time complexity using random generators rand 

V 2 1 

and halton, random sequence sources pi and phi for the integrals 7, = jc dx = — and 

J 3 

o J 

10000 ♦ 

I 2 = [ ^-dx » 1.5707963267949 

io'° * 


Integral with True 
Value 

RSS 

n 

MCI value 

Relative Error 

Time(sec) 

/, « 0.3333333333 


1000 

0.3313420241 

0.0059739277 

0.421 



10000 

0.3331328144 

0.0006015569 

3.813 



100000 

0.3333133046 

0.0000600861 

135.948 



1000 

0.3313370970 

0.0059887089 

0.656 



10000 

0.3331333680 

0.0005998961 

3.891 



100000 

0.3333133335 

0.0000599996 

132.504 


pi 

1000 

0.3313377003 

0.0059868991 

0.438 



10000 

0.3331334104 

0.0005997687 

4.034 



100000 

0.3333133374 

0.0000599877 

752.418 


phi 

1000 

0.3313357881 

0.0059926357 

0.359 



10000 

0.3331332186 

0.0006003443 

3.970 



100000 

0.3333133182 

0.0000600453 

143.128 

/ 2 = 1.5707963268 

rand 

1000 

1.6210367932 

-0.0319840743 

1.494 



10000 

1.5737041052 

-0.0018511492 

14.089 



100000 


0.0001759818 

235.874 j 


halton 

1000 

1.5741008831 

-0.0021037459 

2.415 



10000 

1 .5708375869 

-0.000262670 

14.147 



100000 

1.5708984346 

- 0.0000650039 

239.821 \ 


Pi 

1000 

1.558832916 

0.0075840738 

1.483 



10000 

1.5704023188 

0.0002508333 

14.440 



100000 

1.5708585727 

-0.0000396270 

235.542 


phi 

1000 

1.5905901066 

-0.0126011116 

1.498 



10000 

1.5723190078 

- 0.0009693688 

14.119 



100000 

1.5710505371 

-0.0001618353 

244.999 












































































To conserve space we include in Table 2 Monte Carlo integration values for the following single 
integrals I 3 , U, h, h, h, h, g(0.5), g(l), and g(2) along with relative error in each of rand, halton, pi, 
and phi for the number of subintervals n = 10000 only. 


Table 2 Monte Carlo integration (MCI) along with time complexity using random generators rand 


10000 


r s in x tc 

and Halton, random sequence sources pi and phi for the integrals / 3 = f ( ) 2 dx « — , 

io-'° * 2 

10000 J 1 J 10000 ~ 

7 4 = f -dx*-,I s = f xsin-^r* 0.3785283345, I 6 = 

4 J 1 + x 2 2’ 5 J , 0 x 6 |x 2 +4 4e 6 

10000 ■ 1 

I 7 = f ^ = flnxrft « -1 , / 9 = g(0.5), / I0 = g(l), /„ = g(2) 

« l + x 2e J„ 


Integral with True 
Value 

RSS 

n 

MCI value 

Relative Error 

Time(sec) 

/ 3 = 1.5707963267949 

rand 

10000 

1.5659602830 

0.0030787211 

13.704 


halton 

10000 

1.5709818944 

-0.0001181360 

14.141 


Pi 

10000 

1.5703525786 

0.0002824989 

13.954 


phi 

10000 

1.5722664771 

-0.0009359268 

13.875 

I 4 =l. 5707963267949 

rand 

10000 

1.5763960089 

-0.0035648684 

4.797 


halton 

10000 

1.5707451016 

0.0000326110 

4.719 


Pi 

10000 

1.5701990990 

0.0003802070 

4.828 


phi 

10000 

1.5721051426 

-0.0008332180 

4.562 

I 5 ~0. 3785283345 

rand 

10000 

0.3783614909 

0.0004407690 

10.968 

(true value is 

halton 

10000 

0.3783617259 

0.0004401483 

11.078 

unknown; it is found 

Pi 

10000 

0.3783617763 

0.0004400152 

10.845 

with n=l million) 

phi 

10000 

0.3783616031 

0.0004404728 

11.813 

l 6 =0.0038936 1481 

rand 

10000 

0.0035024365 

-0.7990667613 

15.953 

(sensitive) 

halton 

10000 

0.0019306096 

0.0083202139 

15.935 


Pi 

10000 

0.0019617511 

-0.0076760232 

15.888 


phi 

10000 

0.0021763239 

-0.1178937972 

15.827 

/;= 0.577863674895461 

rand 

10000 

0.5777126672 

0.0002613207 

14.046 


halton 

10000 

0.5778474094 

0.0000281477 

14.187 


Pi 


0.5779883744 

-0.0002157940 

13.812 


phi 

10000 

0.5780012634 

-0.0002380986 

13.687 

h= -0.999999997597415 

rand 


-0.9999905104 

0.0000094872 

12.798 


halton 


- 0.9999998752 

0.0000001224 

12.829 


Pi 


-0.9999996467 

0.0000003509 

12.782 


phi 

F/'W'f'M 

-1.0000019957 

-0.0000019981 

13.016 

I 9 = 1.5707963267949 

rand 

10000 

1.5725361655 

-0.0011076157 

21954 


halton 

10000 

1.5707490321 


21.938 


pi 

10000 

1.5704613149 


22.078 


phi 

10000 

1.5723761285 

-0.0010057330 

21.734 

1,0 = 0.785398163397448 

rand 

10000 

0.7897715914 

-0.0055684215 

21.969 


halton 

10000 

0.7856054911 

-0.0002639778 

21.563 


Pi 

10000 

0.7850052219 

0.0005003087 

21.938 


phi 

10000 

0.7869130037 

-0.0019287546 

21.672 




0.0002371015 

wnmmasm 

22.625 




m — 

Inf 

22.485 


pi 

10000 

-0.0005259715 

Inf 

23.016 


phi 

10000 

0.0013499859 

JsL 

22.563 






























































































We now extract the relative errors and represent these errors in Table 3 for all the eleven typical 
integrals for the sake of ranking the pseudo-random Matlab generator rand, the random sequence 
sources pi and phi, and the quasi-random generator halton in a statistical sense. 

Table 3 (an extraction of Tables 1 and 2) Ranking of pseudo-random generator rand, random 
sequence sources phi and pi, and quasi-random generator halton in terms of relative error for the 
integrals considered here. 


Integral 

rand 

phi 

Pi 

halton 

Best 


rel . error 

rel. error 

rel. error 

rel. error 


Ii 

0.0006015569 

0.0006003443 

0.0005997687 

0.0005998961 

Pi 

I 2 

-0.0018511492 

-0.0009693688 

0.0002508333 

-0.000262670 

Pi 

I 3 

0.0030787211 

-0.0009359268 

0.0002824989 

-0.0001181360 

halton 

I 4 

-0.0035648684 

-0.0008332180 

0.0003802070 

0.0000326110 

halton 

Is 

0.0004407690 

0.0004404728 

0.0004400152 

0.0004401483 

Pi 

16 

-0.7990667613 

-0.1178937972 

-0.0076760232 

0.0083202139 

Pi 

I? 

0.0002613207 

-0.0002380986 

-0.0002157940 

0.0000281477 

halton 

Is 

0.0000094872 

-0.000001998 

0.0000003509 

0.0000001224 

halton 

19 

-0.0011076157 

-0.0010057330 

0.0002132752 

0.0000301087 

halton 

1 10 

-0.0055684215 

-0.0019287546 

0.0005003087 

-0.0002639778 

halton 

111 

infinity 

infinity 

infinity 

infinity 

halton 


Since the true value in integral In is zero (Table 2), the concerned relative error defined as [(true 
value - computed value)/true value] will be evidently infinity as shown in both Tables 2 and 3. 
However, in the absolute sense, the best result was produced by the halton random generator (see 
Table 2). 

From the foregoing result (Table 3), it is clear that the performance of pi as an RSS is consistently 
always better than phi (the golden ratio (p ) without any exception. The transcendental number pi 
{n) is even comparable to the quasi-random generator halton in terms of uniformity. A quasi- 
random generator is known to produce more uniformly distributed random numbers (with less 
discrepancy) than those produced by a pseudo-random generator such as the generator rand. As a 
matter of fact, the digits of pi are more uniformly distributed than those of phi as depicted by the test 
statistic x 2 > viz., the value of tsc (Table 4) and are nearly as uniform as those produced by halton. 
When we consider larger number (more than 10000) of digits, the test statistic j 2 (tsc) is expected to 
be better (smaller) for n than that for <p . Of course, digits fewer than 7500 (say, 2500, 5000) may 
not be sufficient to compare uniformity of digits of these two numbers (p and n . 

In the ranking for the purpose of uniform (at least) one dimensional scanning and producing good 
accuracy, we have halton (rank 1), pi (rank 2), phi (rank 3), and rand (rank 4). It can be observed that 
pi has performed significantly better than both phi and rand in all the foregoing typical examples; it 
has even performed better than halton (and hence others) in about 40% of the problems considered 
here. 

Visualization of the integrands The integrands of some of the integrals considered here are 
interesting and can be visualized so as to appreciate the numerical integration using the Monte Carlo 
methods. We present the graph of these integrands along with the respective Matlab commands. 





“esin x 7t 

Integral I 2 : Visualization of the integrand of I 2 , viz., f dx = — = 1.5707963267949 is 

o * 2 

achieved through the Matlab command »fplot('sin(x)/x', [ 1 0 A - 1 0, 100]) as in Fig. 2. 



Fig. 2 Graph of Limits : [1 O' 10 ,100] 

X 


oo . 

f 1 71 

Integral / 4 Visualization of the integrand of / 4 , viz., (a) I — —^dx = — ; True value= 

0 1 + x 2 

1.5707963267949 is performed by the Matlab command »fplot('l/(l+x A 2)', [0,100]) as in Fig. 3. 



Fig. 3 Graph of -, 

1 + x 


Limits : [0,100] 


Integral / 4 Visualization of the integrand of / 4 , viz., 
the command »fplot('l/(l+x A 2)', [0,5]) as in Fig. 4. 


(b) 


n 

~2 


near origin is done using 





Integrand of / 5 The graph of the integrand (a) jtsin(— ), Limits : [0,1] is produced by the command 

X 

»fplot('x*sin(l/x)', [0,1]) as in Fig. 5. 

i, , . . . , . . . . 1 


0.8 







Fig. 6 Graph of xsin(— ), Limits : [0,0.1] 
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Integrand of I 6 The integrand (a) 


cos 3x 

' 2 +4' 


Limits : [-<»,oo] ; (Integration value = 


n 


Even 

x~ + 4 2e° 

integrand; True numerical value= 0.00389361481414237 ) has the graphical representation as in 
Fig. 7 produced by the command fplot('cos(3*x)/(x A 2+4)', [-100,100]) is 



cos 3x 

Fig. 7 Graph of — , Limits : [-100,100] 

x +4 

cos 3x 

Integrand of / 6 The graph of the integrand (b) — , Limits : [-co,oo] expanded near the origin 

x +4 

and produced by the command »fplot('cos(3*x)/(x A 2+4)', [-10,10]) is shown in Fig. 8. 



cos 

Fig. 8 Graph of — , Limits : [-10,10] 

x +4 



Integral / 7 The integrand of the integral 7, , viz., (a) 

value= 1.15572734979092) is visualized as 
»fplot('x*sin( 1 /x)/( 1 +x A 2)', [-100, 1 00]) 


f — S ' n * dx = — , (Even integrand; True 
i \ + x 2 e 

in Fig. 9 using the command 



x sin x 

Fig. 9 Graph of ~^— r . Limits : [-100,100] 
l + x 


Integral / 7 The integrand of the integral / 7 , i.e., (b) [ - — f dx = — , (Even integrand) in 

' l + jc e 

expanded form near the origin when plotted by the command »fplot('x*sin(l/x)/(l+x A 2)', [- 
0.1,0.!]) would appear as in Fig. 10. 



x sin x 

Fig. 10 Graph of Limits : [-0. 1,0.1] 

1 + JC 



Integral / 8 The integrand of / 8 ,viz., [lnxc/c = xlnx-x = lnl-l-0.011n0.01 + 0.01 = 

0.01 

-0.9439482981401 19 has the graph produced by the command »fplot('log(x)’, [0.01, 1]) as in Fig. 
11 . 



Fig. 1 1 Graph of In x, Limits : [0.01,1] 


Integrals I 9 , 7 10 , 


and/,, The integrands of the integrals g(co) = 


r cos arc sin x 


dx = 


7t !2 , 0 < co < 1 
• 77 1 4 , CO = \ 

0 , w > 1 


Have the graphs shown in Figs. 12, 13, 14. The command »fplot('cos(.5*x)*sin(x)/x', [ 1 0 A - 10,1 00]) 
for co - 0.5 for I 9 depicts the graph as in Fig. 12. 



Fig. 12 Graph of cos ^^ .5xsinx , . j-jq 10 ,100] 
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The command »fplot('cos(x) !|, sin(x)/x', [10^-10,100]) for co = 1 for /, 0 produces the graph as in 
Fig. 13. 





1 



^ ~ cosxsinx . r^-io 

Fig. 13 Graph of , Limits :[10 ,100] 

x 

The command »fplot('2*cos(x)*sin(x)/x , 5 [10 A -10,100]) for a> = 2 for I u generates the graph 
shown in Fig. 14. 



^ COS2xsinX r . . r m-10 mm 

Fig. 14 Graph of , Limits : [10 ,100] 
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Table 4 Distribution of digits of cp and n with test statistic (tsc), critical value of / 2 (cvc= 
14.6840) at 10% significance level and expected frequency (ef) of each digit 



tp (number of 

leftmost digits) 

n (number of 

leftmost digits) 


7500 

10000 

7500 


0 

758 

1020 

701 

968 

1 

795 

1062 

782 

1026 

2 

775 

994 

741 

1021 

3 

784 

1039 

744 

975 

4 

736 

976 

761 

1013 

5 

741 

988 

784 

1046 

6 

672 

918 

762 

1021 

7 

760 

1025 

745 

969 

8 

733 

987 

719 

947 

9 

746 

991 

761 



tsc=14. 1813 

tsc=14.12 

tsc=8.0933 

tsc=9.4580 


ef=750 

ef=1000 

ef=750 

ef=1000 


4. Conclusions 

Study of transcendental and algebraic numbers for random sequence sources for randomized 
algorithms Our numerical experiment shows that pi (x) has always scored over phi (the golden ratio) 
in one dimensional Monte Carlo integration. This implies that other transcendental numbers as well as 
algebraic numbers could also be investigated to determine the uniformity/discrepancy of distribution 
of random numbers generated out of them. This has not only academic interest but also possibly 
significant real world applications. However, while quasi-random sequences are better suited for 
uniform scanning of a space such as that required in Monte Carlo integration, these need not fare 
better than the pseudo-random sequences for other kinds of problems such as the traveling salesman 
problem (TSP) [2], These are because pseudo-random numbers are a bit too random unlike quasi- 
random numbers. Quasi-random numbers, on the other hand, are more uniformly distributed (with less 
discrepancy) than pseudo-random numbers. 

Ready generation of random numbers versus storage-and-retrieval of random numbers Our time 
complexity experiment in Monte Carlo integration demonstrates that unless we have a more efficient 
storage and retrieval system, there is no appreciable advantage of storage and retrieval over the instant 
generation of random numbers. Although storing all the random numbers before the execution of an 
algorithm needs considerable storage space, it may be considered as not an important issue since (i) we 
have sufficient disk space available to us now and (ii) the problem space may be appropriately divided 
into subspaces so that we can always use fixed number (rather small set of) of random numbers for 
each subspace. However, we have experienced significant improvement in time complexity by pre- 
storage and retrieval of random numbers over certain random number generators while solving the 
TSP. 

Use of Pi and Phi in Multiple integration Does pi fare better than phi for double, triple, and higher 
dimensional integrations? Although it might appear that pi should fare better than phi in higher 
dimensional integration as it has done for single integration, we are yet to explore this aspect. This 


















involves significant research mainly on the method of using the random sequence as well as the 
technique of implementation. A detailed study on this will appear elsewhere. 

Connection of pi with Monte Carlo method in history Buffon posed and solved the following problem 
in 1777. Suppose a needle of length L is thrown randomly onto a horizontal plane marked with equi- 
spaced parallel straight lines having a distance d > L between any two consecutive lines. If the 
distance of the center of the needle from the nearest line is x and the angle that its orientation makes 
with the line (or for that matter with any other line) is a , then the needle will intersect a line if and 
only if x < 0.5Z,sina . What is the probability P(x < 0.5Z,sina)that the needle will intersect a line 
[39]? Buffon derived the probability of intersection as 

P = L [sin ad a) !{nd) = — => ;r = — . 

0 J nd dP 

This was an entirely new method of computing n in 1777. Here to calculate the probability P , the 
needle needs to be thrown onto the ruled plane for a large number ( H ) of times and the number (h) of 
times it intersects a line. The probability P = h/H is then computed.. However, during those days 
when today’s modem computer was not even possibly imagined, determining n even to two decimal 
places in this way was not at all attractive; it had definitely an academic interest though. Today with 
the advent of super-high speed computers (over one billion floating point operations per second), such 
a calculation of probability to an acceptable accuracy is not unrealistic. Similar historical relation 
between the golden ratio and the Monte Carlo method is not known to us. 

Deterministic mathematicallydirec methods versus Monte Carlo methods for integration Can the 
deterministic methods such as the closed quadrature formula Trapezoidal rule/Simpson’s 1/3 rule and 
the Gauss-Legendre open quadrature always excel the Monte Carlo methods? The answer is usually 
‘yes’ from both accuracy and computational complexity point of view. The MC procedure is 
sometimes preferred over a numerical quadrature when the integrand is violently fluctuating and/or 
involves long trigonometric/special functions. However, our focus is on the comparison of pi and phi 
— which one is a better source for a random sequence for MC integration/for a uniform scanning of 
given a domain. 

Complexity Issues For most real world applications where the mathematical model consists of a single 
integral, computational/time complexity is often not an important issue since the time of computation 
for four significant digit (0.005%) accuracy is often of the order of seconds. It may be noted that the 
present day (2007) computer can execute over a billion floating operations per second and both Monte 
Carlo and quadrature methods are polynomial-time. From Tables 2 and 3, it appears that there are no 
significant differences in time complexities among the Monte Carlo integrations using any of the four 
RSS’s, viz., rand, halton, pi, and phi. However, for multiple integration, complexity could become an 
issue. We will explore this in our future work. 

Parallel implementation For the MC integration where we divide the region of integration into sub- 
regions, parallel implementation is straightforward since integration for each sub-region can be carried 
out independently and parallely/simultaneously using available processors in a multiprocessor system. 
The outputs are then added in a parallel mode to get the required integration value. 


Appendix 



The following Matlab program mclintegrationrhpiphi stored in the file named 
“mclintegrationrhpiphi” is the Monte Carlo one dimensional integration for all the integrals 
presented in section 3. This program uses another Matlab program halton stored in a separate file. It 
computes the value of an integral for varied number of subintervals with a fixed number N (=5000) 
of random numbers used in each subintervals. It includes the four random number sources, viz., 
rand, halton, pi, and phi and is self-explanatory. Observe that there are four appropriate lines of code 
from which the percentage symbol “%” need to be removed for the execution of the program. To 
conserve space we reproduce, in the Matlab program mclintegrationrhpiphi, tc denoted by piO as 
well as (p denoted phi partially. We have removed the decimal point immediately after the first 
digits 3 and 1 in piO and phi as this is of no concern to us in this random sequence source (SSR) 
context. However, in the actual Matlab program piO and phi must be completely included. We have 
demonstrated later how all the 50000 digits used here for each of piO and phi could be generated 
using the Matlab command v pa and the Matlab program insertblanks. 


^mclintegrationrhpiphi 
function [true, I, comp_accuracy ] =mclintegrationrhpiphi ( fun) ; 

%50000 digits of pi is in the following location piO; each row has 10 elements. 

%Middle two elements are not included so as not to exceed the line length. 

pi0= [3141592653 5897932384 6264338327 9502884197 . . . 4592307816 4062862089 9862803482 5342117067 


9821480865 1328230664 7093844609 5505822317 . . 
6442881097 5665933446 1284756482 3378678316 . . 
3724587006 6063155881 7488152092 0962829254 . . 


8788053304 2714630119 4158989632 8792678327 . . 
6612345888 7310288351 5626446023 6719966445 . . 
9875235781 2331342279 8665035227 2536717123 . . 
8158405229 5369374997 1066559489 4459246286 . . 


. 2841027019 3852110555 9644622948 9549303819 
. 2346034861 0454326648 2133936072 6024914127 
. 0011330530 5488204665 2138414695 1941511609 


. 0466587100 5000832851 7731177648 9735230926 

. 1511493409 3934475007 3025855814 7561908813 

. 6007956982 7626392344 1071465848 9578024140 

. 5339439142 1112718106 9105229002 4657423604]; 


phi=[ 1618033988 7498948482 0458683436 5638117720 . . . 8622705260 4628189024 4970720720 4189391137 
4847540880 7538689175 2126633862 2235369317 . . . 3890865959 3958290563 8322661319 9282902678 

8067520876 6892501711 6962070322 2104321626 . . . 4975870122 0340805887 9544547492 4618569536 

4864449241 0443207713 4494704956 5846788509 . . . 6478091588 4607499887 1240076521 7057517978 


9670572589 3974080648 3175935357 8603833221 
0043613420 1605741421 0103382425 7773560259 
9130860163 5581160937 2610924981 1785171527 
5411334656 3608193864 3097371691 1136979496 

t0=clock; 

^**+********++**********+****+***+****************** 

%FIRST CHANGE: Supply values for alpha, beta, n 

%alpha and beta are lower and upper limits of integration. 

%alpha=0 . 000005; beta= . 01 ; n=1000; a(l)=alpha; 

alpha=10 A -10; beta=10000 ; n=10000; a(l)=alpha; %Interval (beta-alpha) is divided into n subintervals 
%alpha=10 A -10; beta=l; n=10000; a(l)=alpha; %Interval (beta-alpha) is divided into n subintervals. 
%alpha=0; beta=l ; n=10000 ; a(l)=alpha; %Interval (beta-alpha) is divided into 
%n subintervals. This is for integral x A 2 from 0 to 1 . 

%alpha~0 . 000005;beta=l; n~100000; a(l)=alpha; %Interval (beta-alpha) is divided into n subintervals. 
%alpha=0;beta=10000; n=10000; a(l)=alpha; %Interval (beta-alpha) is divided into n subintervals. 
%alpha=0; beta=200; a(l)=alpha; n=10000; %Interval (beta-alpha) is divided into n subintervals. 

for i=2:n-l, a (i) =a ( i — 1 ) + (beta-alpha) /n; end; 

^-a-************************************************** 

N=5000;S1=0; 

********* -k ********************* -k **************** * 

%SECOND CHANGE: Remove '%" from one of rand, halton, piO, and phi 

%x=rand(l, 5000); fprintf ( 'Matlab rand used as random number source\n') %Matlab rand 
% 

%x=halton (1, N) ; fprintf (' halton used as random number source\n') %halton 
% 

%reshape (piO, 1, 5000) ; x=10 A -10*pi0; ; fprintf ( 'pi used as random number source\n') %pi 
% 


. . . 6229114328 2590723291 9935637853 2927911544 

. . . 2438503757 7755389867 5571298529 3795171358 

. . . 9680890339 9517122420 3675938612 4862023009 

. . . 4373853556 2868657173 5845963512 3774246170]; 



reshape (phi, 1, 5000) ; x=10 A -10*phi; fprintf ('phi used as random number source\n') %phi 
^*************************************************** 

fprintf ( 1 # of subintervals n # of random nos. N TrueValue IntegrationValue RelativeError 

Time_in_sec\n ' ) 

for i=2:n-l, S=0;for j-l:N, 
t=x (j); 

^* ******** ************ ******************************************* ************ 
p=a (i) -a (i-1) ;q=a (i-1) ; xl=p*t+q; 

%THIRD CHANGE: Write fun = p* (integrand in terms of xl) 

%fun=p*xl A 2; %Integrand is x A 2 
%fun=p*sin (xl) /xl; %Integrand is sin x/x 
%fun=p* (sin (xl) /xl) A 2; %Integrand is (sin x /x) A 2 
%fun=p* (1/ (l+xl A 2) ) ; %Integrand is l/(l+x A 2) 

%fun=p*xl*sin (1/xl) ; % Integrand is x sin(l/x) 

%fun=p*cos (3*xl) / (xl A 2+4) ; %Integrand is cos 3x /(x A 2+4) 

%fun=p*xl*sin (xl) / (l+xl A 2) ; %Integrand is x sin x/ (l+x A 2) 

%fun=p*log (xl) ; %Integrand is In x 

%fun=p*cos (0 . 5*xl) *sin (xl) /xl; %Integrand is (cos 0.5x sin x) /x 
%fun=p*cos (xl) *sin (xl) /xl; %Integrand is (cos x sin x)/x 
fun=p*cos (2*xl) *sin (xl) /xl; %Integrand is (cos 2x sin x)/x 
^************* ****************************************************** ********* 

S=S+fun; 

%S=S+ (a (i) -a (i-1) ) * ( (a (i) -a (i-1) ) *t+a (i-1) ) * sin (1/ ( (a (i) -a (i-1) ) *t+a (i-1) ) ) ; 
end; 

S=S/N; S1=S1+S; I-Sl ; 
end; 

%************************************************************************** 

% FOURTH CHANGE: Write true (True integration value) 

%true = 1/3; %True value of integral x A 2 from 0 to 1 

%true= 0.3785283345; % Probably true value of integral x sin (1/x) from 10 A -10 to 1 
%true=pi/2; %True value of integral l/(l+x A 2)or sin x/x or (sin x/x) A 2 from 0 to inf. 

%true=-. 943948298140119; %True value of integral In x from .01 to 1 
%true=pi/ (4*exp (6) ) ; %True value of integral cos 3x/(x A 2+4) from 0 to inf. 

%true=pi/ (2*exp (1) ) ; %True value of integral x sin x/(l+x A 2) from 0 to inf. 

%true=l*log ( 1 ) -l-10 A -10*log (10 A -10) +10 A -10; %True value of integral In x from 10 A -10 to 1. 

%true=pi/2; %True value of the integral cos(0.5x)sin x/x from 10 A -10 to inf. 

%true=pi/4; %True value of integral cos x sin x/x from 0 to inf. 
true=0; %True value of integral cos 2x sin x/x from 0 to inf. 

% 

comp_accuracy =(true - I) /true; format long g; 

fprintf (’ %16i %20i %17.10f %17.10f %17.10f %13.6f \n ' , i+1, N, true, I, comp_accuracy, etime (clock, 
t0)); 

The Matlab program “halton” presented below needs to be saved in the Matlab file called halton — a 
file different from the file “mclintegrationrhpiphi”. 

function r = halton (p,q) % where p = 1 and q = n = # of RNs 
dim = 2; persistent seed seed2 base 

if isempty (seed2) , prm_numbers=primes (300) ; base=prm_numbers (dim) ; seed2=base A 4-l ; 
end 

if nargin < 1, p = 1; end, if nargin <2, q - 1; end 

r = zeros (p,q); seed = seed2; 

for k = 1 : p*q, x = 0.0; base_inv = 1.0/base; 

while ( any { seed ~= 0 ) ) 

digit = mod ( seed, base ) ; x = x + digit * base_inv; 
base_inv = base_inv / base; seed = floor ( seed / base ) ; 
end 

r ( k) = x; seed2 = seed2 + 1; seed = seed2; 
end 


Inserting blanks in a string of blankless digits/characters We have used the Matlab program 
insertblanks to convert a .txt file containing consecutive characters/digits without any blank 
anywhere. Consider, for instance, the first 50000 digits of pi or those of phi produced by the Matlab 
variable precision arithmetic command vpa described later in this appendix. The output, viz., the 
string of 50000 digits will be obtained/shown on the Matlab command window in a single line 
without any blank in the string. The insertblanks program takes this string in the form of a .txt file as 
its input and then outputs the string in a matrix with desired number of columns and with desired 
blanks after each fixed-sized block of, say, 10 digits. Insertion of such blanks is necessary to make 



use of the consecutive blocks as random numbers. The program which is self-explanatory is as 
follows. 


%Inserting blanks in .txt file: Blanks are inserted among fixed blocks of (alphanumeric) characters. 

function f = insertblanks() 

clc; 

% show .txt files 
system('dir *.txf); 

% systemfdir •.asv'); % system('dir *.**); 
disp(sprintf( , \n , )); 

inputfile=input('Enter input file name: Vs*); %lnput file has extension .txt 
disp(sprintf(*\n\n')); 

outputfile=input('Enter output file name: ','s'); 
disp(sprintf(’\n')); 

columns=input('Enter no. of columns for display matrix: '); 
dispfsprintfOnVnVn')); 

elementsize=input('Enter no. of digits in each element: ’); 

fldl = fopen(inputfile); 
fid2 = fopenjoutputfile/w'); 

St = fscanfifidl^/os’); loopCount =length(St)/ elementsize; 

S=length(St); initial=1; final=elementsize; Columns = 0; Line = 

for i=1:elementsize:S 

Temp=St(1,initial:final); initial = initial+elementsize; final = final+elementsize; 
if final > S % last number is not K-digits 
final = S; 
end 


Columns = Columns +1; Line = strcat(Line, sprintf(' %s', Temp)); 
if Columns == columns 

disp(sprintf('%s\n',Line)); fprintf(fid2,'%s\n , ,Line); 

Columns = 0; % Reset column counter 
Line-'; 
end 
end 

disp(sprintf('%s\n , ,Line)); % Display left-over 
fprintf(fid2, , %s\n',Line); % Store left-over in file 

fclose(fidl); 

fctose(fid2); 

disp(sprintff\nlnput File Name: %s',inputfile)); 
disp(sprintf('\nOutput File Name: %s',outputfile)); 

disptsprintfOnNumber of coulumns of display matrix: %d', columns)); 
dispjsprintfOnN umber of digits in an element: %d', elementsize)); 

To generate, say, 50000 digits of golden ratio as well as those of pi, one may use the Matlab 
variable precision arithmetic command v pa as follows. 

»vpa((sqrt(5)+l)/2,50000), » vpa(2*asin(l), 50000) 


These will produce 50000 digits each of golden ratio and pi in the Matlab command window. The 
50000 digits will be displayed in one line without any blank anywhere. The blanks are then required 
to be appropriately inserted between fixed blocks of digits. Inserting blanks is achieved by the 
Matlab program insertblanks. Each block will be used as a random number. It is often mapped onto 
an interval, say, [0, 1] in most applications. 

Creating digits and inserting blanks If we, for instance, want to insert a blank after each block of 10 
digits of 50000-digit n , then the procedure could be written as 

5.1 Execute in Matlab command window the command »pi50000=vpa(‘2*asin(l)’, 50000) 

5.2 Copy from the Matlab command window the 50000 digits by clicking “Edit” and then “Copy”. 



S. 3 Open a new M-file and paste these 50000 digits. Remove the decimal point 

S. 4 Click on “File” situated at the left top comer. Go to “Save as”, click and save in pi50000.txt. 

S.5 Now go to Matlab command window and execute the command »insertblanks. At the 
following four prompts provide the requisite file names, number of columns in each line, and 
number of digit in each element of a column. If we want 10 columns in each line and 10 digits in 
each element in a column, and the output file name as pi50000outcoll0blockl0, then the required 
inputs will be provided as follows. 

Prompt Input 

Enter input file name: pi50000.txt 

Enter output file name: pi50000outcoll0blockl0 

Enter number of columns for display matrix: 1 0 

Enter number of digits in each element: 1 0 

Thus we obtain the required output which must be plugged in the program mclintegrationrhpiphi 
before we use pi as an RSS. 
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